subroutine mc_moves(actual_energy)
  use rbf_type
  use variable , only : nn,nip,nodes,q,iseed,mov_attempt,pp_bulk,a
 use random_number2
  use ranftest
  implicit none
  integer imc,nclus,i
  real(dp), dimension(8) :: actual_energy
  real(dp)    , dimension(5)  :: nrandom
  nclus              =  int(dble(NM)*ranf(iseed)) + 1
  mov_attempt(nclus) =  mov_attempt(nclus) + 1
  

450  call RANDOM_NUMBER1( nrandom )

  if ( nrandom(1) >= pp_bulk  ) then
     ! Uniform samplig on the surface nodes
     imc = int( nrandom(2)*dble(nn) ) + 1  ; if ( imc > nn ) imc = nn
!     call surf_moves_mc (a(:,imc),imc,nclus,actual_energy)
  else
     imc = nn + int( nrandom(3)*dble(nip) ) + 1; if ( imc > nodes ) imc = nodes
     call bulk_moves_mc (a(:,imc),imc,nclus,actual_energy)
  end if


!  call surf_moves_mc (q_new)
!  call surface_nano(a_new)
end subroutine mc_moves

subroutine bulk_moves_mc(a_imc,imc,nclus,actual)

  use rbf_type
  use variable , only :dmc,iseed,pi,UldG,vi,invtemp,n_mov_accept,acepta,iseed,da_local,a,wall,ai
  use variable , only : dddx,dddy,dddz,ix1,ddx,ddy,ddz,ddx2,ddy2,ddz2,nn

  use random_number2
  use ranftest
!  use ifport
  implicit none
  integer :: imc,nclus,i,j,k,typ,index,ii,imc_tmp
  real (dp) ,dimension(5) :: a_imc,a2_new,a_old,a_new
  real (dp) ,dimension(3) :: vector,gradimc
  real(dp),  dimension(8) :: energy,actual
  integer ,dimension (6) :: ii_node

 
  real(dp) tr2, tr3, aux1, aux2,aux3, u, v, theta,sqrtu,twopi,valor,increm
  real(dp) :: FLdG, FLdG2, FLdG3,FLdGG, FLdGG2, FLdGG3,free_en_L1_i,gradx,grady,gradz
  real(dp) :: du,grad,ener_el_old,ener_el_new
  logical(lgt)  ::  flag_eigen,flag
  real (dp) :: clock_start,clock_finish,c_ini,c_fin
  real (dp) ::  axn1,axp1,ayn1,ayp1,azn1,azp1,axn,ayn,azn,axp,ayp,azp,energy_surf


  aux1 = 1.d0/sqrt(6.d0); aux2 = sqrt(1.5d0) ;  twopi =2.d0*pi ; energy = 0.d0
  index = nclus
  typ = nclus
  a_old = a_imc
  a_new = a_imc


  increm = dmc(nclus)*( ranf(iseed) - 0.5d0 )
  a_new (nclus) = a_imc (nclus)     + increm
  
  ! check_bound
  ! Boundness for the trace --> Included Correct limits from Orlando's arguments


  tr2    =sum(a_new*a_new) ; if ( tr2 > 2.d0/3.d0 )  return ! flag = .false.

  aux1 = 1.d0/(6.d0*sqrt(2.d0)); aux2 = 2.d0*sqrt(3.d0) ; aux3 = 3.d0*sqrt(3.d0) ; a2_new = a_new**2
             
  tr3 =  aux1*(aux2*a_new(1)**3+aux3*a_new(1)*(-2.d0*a2_new(2)+a2_new(3)-2.d0*a2_new(4)+a2_new(5)) &
       + 9.d0*(a2_new(3)*a_new(4)+2.d0*a_new(2)*a_new(3)*a_new(5)-a_new(4)*a2_new(5)))

  u = tr2/6.d0; v = -tr3/6.d0; theta = v/sqrt( u**3 ) ; sqrtu = -2.d0*sqrt(u)
  
  if ( theta > 1.d0 ) theta = 1.d0
  if ( theta < -1.d0 ) theta = -1.d0
  theta = acos( theta ) 

  vector(1) = sqrtu*cos( ( theta + twopi )/3.d0 ) * 1.5d0
  vector(2) = sqrtu*cos( ( theta - twopi )/3.d0 ) + vector(1)/3.d0
  
  ! Boundness for the tensor and its eigenvalues
  if ( vector(1) < -0.5d0 .or. vector(1) > 1.d0 )  return !  flag = .false.

  valor = ( 1.d0 - vector(1) )/3.d0

  if ( vector(2) < - valor .or. vector(2) > valor ) return ! flag = .false.

  valor = 2.d0*( 3.d0*vector(2)*vector(2)+ vector(1)*vector(1))/3.d0

 

  if ( abs( tr2 - valor ) > 1.d-5 ) then
     print *,' Trace Q2 is not the correct value '
     print *, tr2, valor, abs( tr2 - valor )
!     stop
!     flag = .false.
     return
  end if


!**** Landau de Gennes  ******************


  FLdG =(1.d0/6.d0)* (a_old(typ)-a_new(typ)) *(a_old(typ)+a_new(typ) )*(-3.d0 + ULdG)

  SELECT CASE (typ)
     case (1)
        FLdG2 =((a_old(1) - a_new(1)) * (2.d0* (a_old(1)**2 + a_old(1)*a_new(1) + a_new(1)**2) +  &
             3.d0 * (-2.d0 * a_old(2)**2 +  a_old(3)**2 - 2.d0 * a_old(4)**2 +  a_old(5)**2 ))* ULdG )/(6.d0*sqrt(6.d0))
        FLdG3 =-(0.25d0)*(a_old(1)-a_new(1))*(a_old(1)+a_new(1))*(a_old(1)**2+a_new(1)**2+2.d0*(a_old(2)**2 + a_old(3)**2 + a_old(4)**2 + a_old(5)**2))*ULdG
     case (2)
        FLdG2 = -(a_old(2)-a_new(2)) * (sqrt(3.d0)*a_old(1)*(a_old(2)+a_new(2))-3.d0*a_old(3)*a_old(5) ) * ULdG /(3.d0*sqrt(2.d0))
        FLdG3 =-(0.25d0)*(a_old(2)-a_new(2))*(a_old(2)+a_new(2))*(a_old(2)**2+a_new(2)**2+2.d0*(a_old(1)**2 + a_old(3)**2 + a_old(4)**2 + a_old(5)**2))*ULdG
     case (3)
  
        FLdG2 =( (a_old(3)-a_new(3)) * ((a_old(3)+a_new(3))*( sqrt(3.d0)*a_old(1)+3.d0*a_old(4))+6.d0*a_old(2)*a_old(5))  * ULdG) /(6.d0*sqrt(2.d0))
        FLdG3 =-(0.25d0)*(a_old(3)-a_new(3))*(a_old(3)+a_new(3))*(a_old(3)**2+a_new(3)**2+2.d0*(a_old(1)**2 + a_old(2)**2 + a_old(4)**2 + a_old(5)**2))*ULdG
     case (4)
        FLdG2 =(-(a_old(4)-a_new(4))*(-3.d0*a_old(3)**2+2.d0*sqrt(3.d0)*a_old(1)*(a_old(4)+a_new(4)) + 3.d0*a_old(5)**2)*ULdG)/(6.d0*sqrt(2.d0))
        FLdG3 =-(0.25d0)*(a_old(4)-a_new(4))*(a_old(4)+a_new(4))*(a_old(4)**2+a_new(4)**2+2.d0*(a_old(1)**2 + a_old(2)**2 + a_old(3)**2 + a_old(5)**2))*ULdG
     case (5)
        FLdG2 =( (a_old(5)-a_new(5))*(6.d0*a_old(2)*a_old(3)+(sqrt(3.d0)*a_old(1)-3.d0*a_old(4)) *(a_old(5)+a_new(5)))*ULdG)/(6.d0*sqrt(2.d0))
        FLdG3 =-(0.25d0)*(a_old(5)-a_new(5))*(a_old(5)+a_new(5))*(a_old(5)**2+a_new(5)**2+2.d0*(a_old(1)**2 + a_old(2)**2 + a_old(3)**2 + a_old(4)**2))*ULdG
 END SELECT

energy (1) = FLdG + FLdG2 + FLdG3

!**** Elastic Free Energy  ******************

ii_node(1:6)=  ix1(1:6,imc)
ener_el_old = 0.d0 ; ener_el_new = 0.d0

axn1 = a(typ, ix1(1,ii_node(1)))
axp1 = a(typ, ix1(2,ii_node(2)))
ayn1 = a(typ, ix1(3,ii_node(3)))
ayp1 = a(typ, ix1(4,ii_node(4)))

if  (ii_node(5).gt.nn)then
   azn1 = a(typ, ix1(5,ii_node(5)))
else
   azn1 = 0.d0
endif

if  (ii_node(6).gt.nn)then
   azp1 = a(typ, ix1(6,ii_node(6)))
else
   azp1 = 0.d0
endif

axn = a(typ,ii_node(1))
ayn = a(typ,ii_node(3))
azn = a(typ,ii_node(5))
axp = a(typ,ii_node(2)) 
ayp = a(typ,ii_node(4))
azp = a(typ,ii_node(6))


 if(.not.wall(imc))then


    energy (2)=    (a_new(nclus) - a_old(nclus))*(-((axn1 + axp1)*dddx**2) -  (ayn1 + ayp1)*dddy**2 - (azn1 + azp1)*dddz**2 & 
         + a_new(nclus)*(dddx**2 + dddy**2 + dddz**2) +   a_old(nclus)*(dddx**2 + dddy**2 + dddz**2))


!    if(ranf(iseed).lt.0.5d0)then

       ener_el_old = ((a_new(nclus) - a_old(nclus))*((3.d0*a_new(nclus) + 3.d0*a_old(nclus) + 2.d0*(-6*azn + azn1 + 2*azp))*ddx2**2*ddy2**2 + ((3.d0*a_new(nclus) & 
            + 3.d0*a_old(nclus) + 2.d0*(-6.d0*ayn + ayn1 + 2.d0*ayp))*ddx2**2 + (3.d0*a_new(nclus) + 3.d0*a_old(nclus) + 2.d0*(-6.d0*axn + axn1 + &
            2*axp))*ddy2**2)*ddz2**2))/(24.d0*ddx2**2*ddy2**2*ddz2**2)


       ener_el_new =((a_new(nclus) - a_old(nclus))*((3.d0*a_new(nclus) + 3.d0*a_old(nclus) + 2.d0*(2.d0*azn - 6.d0*azp + azp1))*ddx2**2*ddy2**2 + ((3.d0*a_new(nclus) &
            + 3.d0*a_old(nclus) + 2.d0*(2.d0*ayn - 6.d0*ayp + ayp1))*ddx2**2 + (3.d0*a_new(nclus) + 3.d0*a_old(nclus) + &
            2.d0*(2.d0*axn - 6.d0*axp + axp1))*ddy2**2)*ddz2**2))/(24.d0*ddx2**2*ddy2**2*ddz2**2)


    energy( 2) = energy (2) + (ener_el_new  +ener_el_old)*0.5d0
    energy_surf = 0.d0

 else


    if (ii_node(5).gt.nn)then  
       energy(2) = ((a_new(nclus) - a_old(nclus))*(2.d0*(a_new(nclus) + a_old(nclus) - axn1 - axp1)*dddx**2 + 2.d0*(a_new(nclus) + a_old(nclus) - ayn1 - ayp1)*dddy**2 &
            + (a_new(nclus) + a_old(nclus) - 2.d0*azn1)*dddz**2))/2.d0

       energy_surf =   0.5d0*(((a(index,ii_node(6)) - a_new(nclus))  /ddz  ) **2- ((a(index,ii_node(6)) - a_old(nclus))  /ddz  ) **2)

    else if (ii_node(6).gt.nn)then 
       energy(2) = ((a_new(nclus) - a_old(nclus))*(2.d0*(a_new(nclus) + a_old(nclus) - axn1 - axp1)*dddx**2 + 2.d0*(a_new(nclus) + a_old(nclus) - ayn1 - ayp1)*dddy**2 &
            + (a_new(nclus) + a_old(nclus) - 2.d0*azp1)*dddz**2))/2.d0

      energy_surf =   0.5d0*(((-a(index,ii_node(5)) + a_new(nclus))  /ddz  ) **2- ((-a(index,ii_node(5)) + a_old(nclus))  /ddz  ) **2)
    endif

    ener_el_old = ((a_new(nclus) - a_old(nclus))*((3.d0*a_new(nclus) + 3.d0*a_old(nclus) + 2.d0*(-6*azn + azn1 + 2*azp))*ddx2**2*ddy2**2 + ((3.d0*a_new(nclus) & 
         + 3.d0*a_old(nclus) + 2.d0*(-6.d0*ayn + ayn1 + 2.d0*ayp))*ddx2**2 + (3.d0*a_new(nclus) + 3.d0*a_old(nclus) + 2.d0*(-6.d0*axn + axn1 + &
         2*axp))*ddy2**2)*ddz2**2))/(24.d0*ddx2**2*ddy2**2*ddz2**2)

    ener_el_new =((a_new(nclus) - a_old(nclus))*((3.d0*a_new(nclus) + 3.d0*a_old(nclus) + 2.d0*(2.d0*azn - 6.d0*azp + azp1))*ddx2**2*ddy2**2 + ((3.d0*a_new(nclus) &
         + 3.d0*a_old(nclus) + 2.d0*(2.d0*ayn - 6.d0*ayp + ayp1))*ddx2**2 + (3.d0*a_new(nclus) + 3.d0*a_old(nclus) + &
         2.d0*(2.d0*axn - 6.d0*axp + axp1))*ddy2**2)*ddz2**2))/(24.d0*ddx2**2*ddy2**2*ddz2**2)


   
    energy( 2) = energy (2) +( ener_el_new + ener_el_old)*0.5d0 

endif


energy (1) = energy (1)* vi
energy (2) = energy (2)* vi + energy_surf *ai
energy (8) = energy (1) + energy(2)       

   du  =energy (8) * invtemp
 !  write(*,*)energy (8) * invtemp,energy (8) , invtemp

  if ( -du .lt.-30.d0   ) then
     a_imc = a_old
  else if ( -du .ge.0.d0   ) then
     n_mov_accept(nclus) = n_mov_accept(nclus) + 1
     acepta = acepta + 1
     actual =  energy  +   actual
     a_imc = a_new
  else if ( ranf(iseed) < exp( -du )  ) then
     n_mov_accept(nclus) = n_mov_accept(nclus) + 1
     acepta = acepta + 1
     actual =  energy  +   actual   
     a_imc = a_new
  else 
     a_imc = a_old
  end if

!  a_imc (nclus) = a_imc (nclus) - increm

end subroutine bulk_moves_mc
!------------------------------------------------------------------
!------------------------------------------------------------------
!------------------------------------------------------------------
!subroutine surf_moves_mc(a_imc,imc,nclus,actual)

!  use rbf_type
!  use variable , only :dmc,iseed,pi,UldG,vi,invtemp,n_mov_accept,acepta,iseed,da_local,a,wall
!  use variable , only : dddx,dddy,dddz,ix1,ddx,ddy,ddz,ddx2,ddy2,ddz2,nn

!  use random_number2
!  use ranftest
!  use ifport
!  implicit none
!  real(dp), dimension(5) :: v1_gal
!  real(dp), dimension(5,3) :: grad, v2_gal
!**** Rapini  ******************
!  i=node     
  
 ! Surface Energy
  ! Normal vector
!  n   = normal(3*i - 2:3*i)
     
  ! Galatola vectors
!  v1_gal(1) = aux2*( 2.d0*n(3)*n(3)/3.d0 - n(1)*n(1)/3.d0 - n(2)*n(2)/3.d0 )
!  v1_gal(2) = sqrt(2.d0)*n(1)*n(2)
!  v1_gal(3) = sqrt(2.d0)*n(1)*n(3)
!  v1_gal(4) = ( n(1)*n(1) - n(2)*n(2) )/sqrt(2.d0)
!  v1_gal(5) = sqrt(2.0d0)*n(2)*n(3)

!  v2_gal = 0.d0
!  v2_gal(1,1) = -n(1)/3.d0   ; v2_gal(1,2) = -n(2)/3.d0 ; v2_gal(1,3) = 2.d0*n(3)/3.d0 ; v2_gal(1,:) = v2_gal(1,:)*aux2
!  v2_gal(2,1) = n(2) ;         v2_gal(2,2) = n(1) ;                                      v2_gal(2,:) = v2_gal(2,:)*sqrt(2.d0)/2.d0
!  v2_gal(3,1) = n(3) ;                                    v2_gal(3,3) = n(1) ;           v2_gal(3,:) = v2_gal(3,:)*sqrt(2.d0)/2.d0
!  v2_gal(4,1) = n(1) ;         v2_gal(4,2) = -n(2) ;                                     v2_gal(4,:) = v2_gal(4,:)/sqrt(2.d0)
!  v2_gal(5,2) = n(3) ;       v2_gal(5,3) = n(2) ;           v2_gal(5,:) = v2_gal(5,:)*sqrt(2.d0)/2.d0

!  rapini = 0.d0; v_gal = 0.d0; tr_gal = 0.d0

!  do j = 1, 5
!     rapini = rapini + ( a(j,i) - ao_s(j,i) )**2
!     tr_gal = tr_gal + a(j,i)*v1_gal(j)
!     v_gal  = v_gal  + a(j,i)*v2_gal(j,:)
!  end do

!  tr_gal = A5 - tr_gal 
!  r        = dot_product( n, v_gal )
!  galatola = 2.d0*dot_product( v_gal, v_gal ) + 2.d0*r*r + 4.d0*r*tr_gal + tr_gal*tr_gal
     
!  energy(6) =  0.5d0*gamma_w1*( 1.d0 - chi(i) )*galatola + 0.5d0*gamma_s*chi(i)*rapini

!end subroutine surf_moves_mc
